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Abstract The continuous-time random walk (CTRW) 
model is useful for alleviating the computational bur¬ 
den of simulating diffusion in actual media. In prin¬ 
ciple, isotropic CTRW only requires knowledge of the 
step-size, Pj, and waiting-time, Pj, distributions of the 
random walk in the medium and it then generates pre¬ 
sumably equivalent walks in free space, which are much 
faster. Here we test the usefulness of CTRW to mod¬ 
elling diffusion of finite-size particles in porous medium 
generated by loose granular packs. This is done by first 
simulating the diffusion process in a model porous medium 
of mean coordination number, which corresponds to 
marginal rigidity (the loosest possible structure), com¬ 
puting the resulting distributions Pi and Pt as functions 
of the particle size, and then using these as input for a 
free space CTRW. The CTRW walks are then compared 
to the ones simulated in the actual media. 

In particular, we study the normal-to-anomalous tran¬ 
sition of the diffusion as a function of increasing par¬ 
ticle size. We find that, given the same Pi and Pt for 
the simulation and the CTRW, the latter predicts in¬ 
correctly the size at which the transition occurs. We 
show that the discrepancy is related to the dependence 
of the effective connectivity of the porous media on the 
diffusing particle size, which is not captured simply by 
these distributions. 

We propose a correcting modification to the CTRW 
model - adding anisotropy - and show that it yields 
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good agreement with the simulated diffusion process. 
We also present a method to obtain Pi and Pt directly 
from the porous sample, without having to simulate an 
actual diffusion process. This extends the use of CTRW, 
with all its advantages, to modelling diffusion processes 
of finite-size particles in such confined geometries. 

Keywords Anomalous diffusion • Granular pore 
space • Continuous-time random walk • Anisotropy 

1 Introduction 

Diffusion plays a key role in a wide range of natural and 
technological processes. A textbook modelling of such 
processes is the consideration of the diffusion of a single 
memory-free particle in a given medium. The nature of 
such a random walk is governed by three probability 
density functions (PDFs): of the step size, Pi{li); of the 
step direction, Pnint); and of the waiting time between 
steps, Pt[ti). These PDFs are, in principle, position de¬ 
pendent, but it is standard practice to derive (or pos¬ 
tulate) them assuming position-independence and that 
Pnifii) is uniform. The diffusion is then modelled as a 
continuous-time random walk (CTRW) in free space. 
Specifically, the CTRW is constructed by adding vec¬ 
tors of uniformly random orientations, whose lengths 
are chosen from Pi, at time intervals chosen from Pt. 
Averaging over sufficiently many such independent pro¬ 
cesses, the dependence of the mean square distance 
(MSD) on time satisfies (x^) = In normal diffusion 
0 = 1 and D is the standard diffusion coefficient. But 
when Pi and/or Pt are very wide, the diffusion might 
become anomalous (a 7 ^ 1). In particular, when Pt has 
a slowly decaying algebraic tail and Pi does not, the 
random walk is sub-diffusive (a < 1) [B H]- Alterna¬ 
tively, if Pi has a slowly decaying algebraic tail and Pt 
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does not, the random walk is super-diffusive (a > 1), 
resembling a Levy flight [3]. Diffusion processes that 
have the same value of a are said to be in the same 
universality class [1]. 

Anomalous diffusion can arise from different sources, 
which can only be identified by going beyond the MSD. 
When single particle tracking is possible, the movement 
can be evaluated by the time-averaged MSD (TAMSD), 
S‘^{t,T) [3 -While the MSD is the ensemble average of 
the squared distance, made during a time interval t, 
over different realisations, the TAMSD, i5^(t,T), is the 
average of the same quantity along a single trajectory 
of length T. Within the model of sub-diffusive CTRW, 
the TAMSD satisfies (5^) ~ t ■ T““^, where the angu¬ 
lar brackets denote a further ensemble average. In con¬ 
trast, the MSD is sub-linear in t, which makes CTRW 
non-ergodic - the time-average and ensemble-average 
differ. In particular, the dependence of the TAMSD on 
T points to the ageing nature of CTRW [5]. 

A key feature of sub-diffusive CTRW is the random¬ 
ness of its TAMSD. Since P* is scale free, the longest 
waiting times each individual trajectory encounters vary 
significantly, as do the amplitudes of the individual 
TAMSDs. To quantify this, we define the amplitude 
scatter, ^ = 5"^/{5'^). For ergodic processes (e.g. a = 1) 
its PDF is P{^) = S{^ — 1) for sufficiently long tra¬ 
jectory times. But within CTRW this PDF broadens 
as a decreases. Defining the ergodicity breaking (EB) 
parameter, EB = — (^)^, it can be derived analyt¬ 

ically for CTRW processes as a monotonically decreas¬ 
ing function of a. 

Another cause for sub-diffusion is walking in a fractal¬ 
like environment B El- Such environment is charac¬ 
terised by a network of narrow passages and dead ends 
at different length scales, which hinder the walk. Unlike 
CTRW, this process is stationary and therefore ergodic. 
The TAMSD, like the MSD, is sub-linear in t, indepen¬ 
dent of T and its EB parameter vanishes. 

Using CTRW to model diffusion in confined geome¬ 
tries, such as porous media formed by either sintered 
or unconsolidated granular materials, is very attractive 
[i 1 [H [U because it alleviates the need to simu¬ 
late directly the dynamics of particles within the pore 
space, reducing significantly the computational burden. 
In addition, it alleviates finite-size errors due to finite 
samples. This practice is based on the common assump¬ 
tion that Pi , Pt and P„ alone control the random walk’s 
universality class. The common procedure is to find first 
the forms of these distributions in a specific medium, 
using either small simulations or analytic derivation un¬ 
der some assumptions, and then use these to carry out 
a many-step CTRW in free space. It is then presumed 


that the CTRW yields the same universality class as 
the diffusion in the confined geometry. 

The first aim of this paper is to demonstrate that 
this does not apply when the size of the diffusing parti¬ 
cle is comparable to throat sizes. We do so by analysing 
trajectories of individual particles diffusing in a porous 
sample and show statistical deviations from CTRW pre¬ 
dictions. We also compare these simulations with an 
equivalent CTRW model. We show that the effective 
change in the medium’s connectivity with varying par¬ 
ticle size affects directly the nature and universality 
class of the diffusion process. We conclude that the 
sub-diffusion is the result of CTRW on a percolation 
cluster. Indeed, a combination of underlying mecha¬ 
nisms, leading to sub-diffusion, has also been observed 
in [Ulllllllig. 

The second aim of the paper is to propose a method 
to correct for the topological effect, which makes it pos¬ 
sible to still use CTRW, with its advantages, to model 
diffusion of any finite size particle in confined geome¬ 
tries. To maximise the range of validity of our results 
(see discussion below), we consider very high porosity 
porous media. These correspond to marginally rigid as¬ 
semblies of frictional particles, whose mean coordina¬ 
tion number is four m- The least confined of these 
are model systems whose each particle has exactly four 
contacts. 

The structure of this paper is the following. In sec¬ 
tion we describe the simulated porous samples. In 
section we describe the diffusion process, and dis¬ 
cuss the effects of particle size. We perform statistical 
analysis of the particle trajectories and show disagree¬ 
ments with some predictions of the CTRW model. In 
section]^ we describe the equivalent CTRW simulations 
and show that they yield different behaviour in spite of 
having the same step-length and waiting-time distribu¬ 
tions. We propose an explanation for this discrepancy. 
In section we propose a modification to the conven¬ 
tional CTRW model to alleviate this problem, making 
it more suitable for modelling diffusion of finite size par¬ 
ticles in confined geometries. We conclude in section 
with a discussion of the results. 

2 The porous sample 

To simulate a three-dimensional porous granular as¬ 
sembly of coordination number four, we first generate 
an open-cell structure, using Surface Evolver as follows 
EZi ng. Initially, N seed points are distributed ran¬ 
domly and uniformly within a cube, and the cube’s 
space is Voronoi'-tessellated to determine the cell as¬ 
sociated with each point. A cell around a point consists 
of the volume of all spatial coordinates closest to it. 
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The resulting cellular structure is then evolved with 
Surface Evolver to minimise the total surface area of 
the cell surfaces. This procedure is used commonly to 
model dry foams and cellular materials whose dynam¬ 
ics are dominated by surface tension. The result is an 
equilibrated foam-like structure, comprising cells, mem¬ 
branes, edges and vertices. A membrane is a surface 
shared by two neighbouring cells, an edge is the line 
where three membranes coincide, and a vertex is the 
point where four edges coincide. 



Fig. 1: Pseudo-grains around the foam vertices. 


3 Diffusion in the porous sample 

We model the diffuser as a sphere of radius r, mea¬ 
sured in units of the average effective throat radius, rg. 
We start by considering particles that are considerably 
smaller than the smallest throat in the structure. The 
particle cannot enter the tetrahedral pseudo-grains, but 
only move from cell c to cell c' through their shared 
throat. The simulation progresses by moving the par¬ 
ticle from one cell, c, to a neighbouring cell, c'. Each 
such an event is a step, lc,c', namely a vector extending 
between the centres of these cells. We define a wait¬ 
ing time, tc, which is the number of time steps spent 
in cell c before a jump occurs. Inside a cell, the parti¬ 
cle is assumed to undergo Brownian motion and is 
proportional to (i) the square of the effective cell ra- 
dius, Rc = (^) j where Vc is the cell volume, and 
(ii) the inverse of the fraction of the cell’s open surface 
through which the particle can pass to neighbouring 

^(t) 

cells. Sc = ■ This is since the particle, on av¬ 

erage, has to make 1/S'c journeys of length Rc until 
it goes through a throat rather than hits a facet of a 
pseudo-grain. The effective area of a throat is the area 
through which a particle of radius r can pass and 

is the total area of throats (facets) that make the 
surface of the cell. is then the sum of the effective 
throat areas accessible for the particle to go through. 
Thus, the waiting time within a cell is 


Next, we construct a tetrahedron around every ver¬ 
tex by connecting the mid-points of the four edges em¬ 
anating from it |19j . Neighbouring tetrahedra are in 
contact in the sense that they share the mid-point of 
an edge. This construction results in a pseudo-granular 
structure of volume fraction (j) = 34%, in which ev¬ 
ery tetrahedron represents a pseudo-grain in contact 
with exactly four others [20] (see hg. [^. Since neigh¬ 
bouring pseudo-grains share the mid-point of the edge 
between them, the tetrahedra structure is topologically 
homeomorphic to the original structure. The void space 
surrounded by the pseudo-grains is still cellular, but 
a cell surface now consists of triangular facets - the 
faces of the pseudo-grains surrounding it ~ and throats 
- the skewed polygons remaining of the original cell 
membranes. The membranes over the throats are dis¬ 
regarded, resulting in an open-cell porous structure, in 
which the throats are the openings between neighbour¬ 
ing cells. The pseudo-grains volumes are smaller than 
those of real convex grains, which curve out into the 
cells of this structure. This increases the pore volume 
and forms a limiting case, which establishes the valid¬ 
ity of our results for any porous medium, as will be 
discussed in the concluding section. 



2 d \ 47r y 


1 + 


A, 


(t) 
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where d is the local diffusion coefficient, which, using 
the Stokes-Einstein relation, is inversely proportional 
to the particle radius, d = (ro/rjdg- 

The probability to exit cell c into c', P(c' | c), is pro¬ 
portional to the area of the throat between them. To 
reduce finite-size effects, we wrap the sample around 
with periodic boundary conditions and let the particle 
travel larger distances by re-entering the sample. This 
means that the same cell may occur at different loca¬ 
tions along the random walk. To avoid distorting the 
statistics by using the same cell too many times, we 
stop the process once a cell has occurred at more than 
five different locations. 

We emphasise that we do not simulate the diffusion 
within the cell - each step in our simulation corresponds 
to a transition of the particle from one cell to another, 
and the waiting time associated with the step is calcu¬ 
lated from the cell properties. This process constitutes 
a random walk on a graph, whose nodes are the cell 
centres. After waiting for a period of time tc in cell c. 
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determined by eq. Q, the particle makes a step to a 
neighbour cell, according to P{c' \ c). 

Before continuing, it is instructive to put the prob¬ 
lem in thermodynamic context. The cell can be re¬ 
garded as a potential well of height AE, and the proba¬ 
bility to escape from it is P{t) = We can then 

use Kramer’s escape rate formula, 

^ ^ 27rfcT ^AE/kT 

d^yU” ia)U” (b) 

where k is the Boltzmann constant, T is the tempera¬ 
ture and U” (a) and U” (b) are the second derivative of 
the potential at the bottom and top of the well, respec¬ 
tively. Interpreting tc as the half-life of the particle in 
the cell, tc = Tln2, we can combine eq. 0 and 0 to 
get: 

21 n(2)7rfcT _ 1 ,^3^^ 

^/U”{a)U”{b) ‘2Sc{r) {in) ' ^ ^ 


Assuming that all cells have the same effective potential 
[/, we get: 

AE 

— = Co«t.-i„r + i„^ . (4) 

This equation establishes the height of the effective bar¬ 
rier in terms of the cell volume and the fraction of its 
surface through which the particle can escape. 

Note that the particle’s mean free path within a cell 
is assumed to be well smaller than the cell size, regard¬ 
less of the particle size, and therefore that Knudsen dif¬ 
fusion [mill] need not be considered. However, even if 
this assumption is not borne out, this would only mod¬ 
ify the coefficient d in eq. Q, which is arbitrary anyway 
in our simulations. Also note that the above assump¬ 
tion, P{t) ~ 6“*/'’’, means that typical escape times do 
not deviate much from the mean or half-life time. This 
justifies our choice of taking p as a representative. 









(e) 


(f) 


(g) 


(h) 


Fig. 2: Statistical analysis of diffusion in the porous sample. All averages are over 1000 walks, (a) MSD vs. time for a small 
particle (r = ro/100), for which a ~ 1, indicating normal diffusion. Inset: overall waiting-time distribution of the walks, (b) 
Step-length distribution, calculated directly from the sample (see text), with a normal fit. Pi is nearly constant for 

all particle sizes, as the cells positions are fixed. P; is narrow, and thus does not affect the universality class of the walk, (c, 
d) Waiting-time distributions (left) and MSD vs. time (right) for large particles (r/ro = 35). Both a and /3 decrease 

with increasing r. All power-law relations hold for at least 1.5 orders of magnitude, (e, f) Ensemble-averaged TAMSD, (5'^), 
vs. t (left) and vs. T (right) for the same large particles. For t > 10^ ®, (5^) becomes sub-linear in t, like the MSD, and is 
nearly independent of T. (g) PDF of the amplitude scatter, from 2000 individual walks of three different particle sizes. The 
distribution broadens dramatically with particle size, which indicates the walk is non-ergodic. (h) The EB parameter for the 
same three particle sizes, as well as for r = ro/100, vs. the anomaly parameter, a. The solid line is the expected relation for 
CTRW, which goes through all points. Error bars denote 2fT ranges, and were established by the variance of 10 independent 
measurements. 
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For each walk we calculate the particle’s position 
at time t, relative to the origin, x(t) = J2n=i 1"; where 
N{t) is the number of steps made before time t. We then 
calculate the MSD, (x^), as a function of t, where the 
angular brackets denote average over 1000 walks. Figure 


The linear relation indicates normal diffusion, with a 
diffusion coefficient oi D = (0.65±0.03)(i {d = lOOdg)- 
The inset figure shows a narrowly bounded Pt- 


2a shows {'xr) vs. t for a particle of size r = tq/IOO. 


Large particles We next consider particles of sizes com¬ 
parable to tq. Such particles diffuse differently due to 
two effects. One is delay and trapping inside cells. The 
larger r is the lower the probability of passing through 
any particular throat, since the effective area, which the 
particle can pass through, is smaller. This reduces the 
overall probability to exit a cell, increasing the waiting 
times spent inside cells. As a result, while the waiting 
times are narrowly bounded for a small particle, Pt{ti) 
develops a power-law tail for sufficiently large particles, 
Pt(ti > ~ , with /3 a function of r. The second 

effect is that the topology changes with particle size; 
as it increases, the probability of passing through some 
throats vanishes identically, changing the system’s con¬ 
nectivity for this particle. 

Fig. shows Pt for a few large particles. We see 
that j3 decreases with r, corresponding to longer wait¬ 
ing times. Fig. |2d| shows the MSD vs. time for the same 
particles. We see that beyond a certain particle size 
(r ^ 1.2ro) a starts to decrease - the diffusion becomes 
anomalous. To quantify this relation, we choose a larger 
set of radii and plot a{r) vs. (3{r) (blue dots in fig. [^. 
We see that for smaller particles /3 ^ 2 and a = 1, 
corresponding to normal diffusion with narrow waiting¬ 
time PDFs. As r increases, f3 decreases and the walks 
eventually become sub-diffusive with a < 1. We mea¬ 
sure a transition at = 2.53 ± 0.03. 

A short comment on scaling windows is due - a is 
evaluated along t G (10^-®, 10®) (see fig. 2d I. At much 
longer times, the diffusion is normal for all particle sizes. 
This is merely a consequence of the finite system size - 
there are no cells with longer waiting times than ~ 10®. 

To investigate the diffusion process further, we cal¬ 
culate the TAMSD, 6“^, and its average over 1000 walks, 
((5^). Fig. 2e andshow {S'^{t,T)) vs. the time lag, t, 
and vs. the overall trajectory time, T, respectively. (5^) 
is sub-linear in t, deviating from the linear t-dependence 
in the CTRW model, and it is nearly independent of T. 
This behaviour is the same as for a random walk on 
a fractal. In contrast, the amplitude scatter, follows 
the CTRW prediction: Fig. shows the PDF of ^ for 
three different particle sizes of order tq. The dramatic 



Fig. 3: Two sets of simulations - diffusion in a porous sample 
(blue dots) and CTRW in unconfined space (green triangles). 
Each simulation (dot) represents a particular particle size, 
and is described by the power-law of the waiting-time dis¬ 
tribution, di and the anomaly parameter, a. Small (large) 
particles appear at the top right (bottom left) corner of the 
graph - they experience a narrow (wide) waiting-time dis¬ 
tribution and undergo normal (sub-) diffusion. The blue and 
green lines are fits of the form a = 1— ^ exp{—(/3 —/3t — ^)/t}, 
with = 2.53, r = 0.42) and (/jp^RW) ^ 2.01, 

T = 0.40), respectively. The theoretical CTRW prediction of 
a universality class transition at /3 = 2 is denoted by the 
black dashed line. The red lines mark the range of /3’s that 
corresponds to the percolation threshold of the sample. This 
range matches the universality-class transition for confined 
diffusion. Error bars denote 2 (t ranges, and were established 
by the variance of 20 independent measurements. 


broadening of the PDF as r increases indicates that 
the process is non-ergodic and that there is ageing [3. 
Fig. 1^ shows that the corresponding EB parameters, 
i.e. the variance of these PDFs, follow the theoretical 
expectation from CTRW. 


4 Modelling the diffusion as isotropic CTRW 

Next we show that an attempt to simulate this pro¬ 
cess with straightforward isotropic CTRW fails. This 
may not come as a surprise, as the statistical analysis 
showed some deviations from the traditional CTRW, 
but it is still constructive to describe the CTRW simu¬ 
lation to better understand its adjustments in sectionj^ 
For such a simulation we use the r-dependent PDFs, Pi 
and Pt, that the particle experiences while diffusing in 
the confined structure. Recall that these PDFs refer to 
the transition between cells, rather than the movement 
within a cell. One way to obtain these PDFs is to mea¬ 
sure them empirically during a diffusion process. How- 
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ever, more efficient is to compute them directly from 
the structural statistics of the porous medium, as we 
outline next. 

A derivation of Pi and Pt from structural statistics 
should be made cautiously because a straight-forward 
histogram of the waiting times of all cells in the sam¬ 
ple ignores the inherent correlation between the prob¬ 
ability to visit a cell, P{c), and the waiting time [23] . 
Cells that are difficult to get out of (long waiting times) 
tend to have a lower probability of getting into and are 
therefore visited less frequently. In particular, some cells 
are completely inaccessible for particles above a certain 
size. 

To this end we use the observation that, to a very 
good accuracy in our diffusion process, P{c) is linear 
in the cell’s total effective throat area, A^\ This ob¬ 
servation, which is independent of the cell volume, can 
be seen over 3.5 orders of magnitude in fig. Further¬ 
more, this holds for all particle sizes, both well smaller 
and comparable to tq. This allows us to estimate P{c) 
for a particular particle size by using the effective 
- a direct structural characteristic of the medium. In 
principle, one expects the visiting probability to be cor¬ 
related with the visiting probabilities of neighbouring 
cells, but fig. shows that this effect is negligible. 



Fig. 4: Visiting probability, P(c), vs. the cell’s total effective 
throat area, in a diffusion process with r = r^. 


Collecting the statistics of all possible steps in the 
sample, we get a PDF described well by the Gaussian 
Pi{li) ^ exp{ —(Zi —r)^/2(T^}, with ajl = 0.15±0.02 (see 
fig. 2b). As Pi is narrow, it does not affect the universal¬ 
ity class of the walk. Pi stays narrow, and indeed nearly 
constant, for all particle sizes, as the cells positions are 
fixed. Note that it is also possible to derive Pi analyt¬ 
ically, given the cell volume distribution and nearest 
neighbour volume-volume correlations. This, however, 
is somewhat downstream from the main thrust of this 
paper. 

We are now able to obtain accurate step-length and 
waiting-time PDFs for every particle size, which could 
be used as input into unconstrained CTRW. Fig. 
shows results of CTRW simulations (green triangles) 
using these PDFs. The fit to this curve (solid green line) 
differs from the theoretical prediction [ 21 ], a = (3 — 1, 
due to finite time and size effects. Our measured transi¬ 
tion for CTRW is = 2.01 ± 0.02, in agreement 

with the theoretical prediction. 

A key observation is that the CTRW exhibits a 
normal-to-anomalous transition at a lower values of /3 
than in the actual sample. Since both processes have 
the same step-length and waiting-time distributions, 
but one is performed on a graph and the other in free 
space, the discrepancy must stem from the connectivity 
of the sample, which the CTRW model cannot account 
for. This is supported by the statistical analysis of the 
diffusion in the porous sample (section , that show 
behaviours typical to random walks on a fractal. 

As mentioned above, the size of the diffusing particle 
determines the effective throat sizes, and hence the con¬ 
nectivity of the porous sample. Moreover, above some 
size there is no path percolating between the sample’s 
boundaries. Fig. shows the dependence of the per¬ 
colating accessible volume on r, where a range of radii 
around the percolation threshold is marked by red lines. 
The same range is marked in fig. We see that the uni¬ 
versality class transition in the sample occurs within 
this range. This is more evidence that the connectivity 
plays an important role in determining the universal¬ 
ity class. Specifically, around the percolation threshold 
the incipient cluster assumes a fractal-like structure, 
further inducing sub-diffusion EHTj. We conclude that 
our simulations of diffusion in porous media are best 
described as CTRW on a percolation cluster. 


We can now estimate more accurately Pt-, and, par¬ 
ticularly, the power-law /3, for every particle size, by ma¬ 
nipulating the waiting-time histogram as follows. For 
every bin consisting of the waiting times of cells {ci,..., c„} 
we multiply the bin’s height by P(ci) and normalise 
the histogram into a PDF. This suppresses long waiting 
times in Pt. 


5 Anisotropic CTRW 

To retain the usefulness of the CTRW model, it would 
be desirable to modify it to capture the effect of connec¬ 
tivity. We next propose such a modification, inspired by 
the PDF of the angle between successive steps in the 
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Fig. 5: The percentage of cells belonging to the incipient- 
cluster vs. the particle radius (in units of rp). The percolation 
range is marked in red. The corresponding 0 range is marked 
in hg. 1^ 


porous sample (fig. |^. There is a finite probability to 
make a backward step, i.e. go back through the throat 
the particle entered a cell, Pback = Pe{d = n). For very 
small particles Pback = 0.087 ± 0.001. This is in con¬ 
trast to the conventional CTRW model, where the next 
step direction is uniformly distributed. Moreover, we 
see that Pback > 1/13.7 « 0.073, which is the inverse of 
the average number of throats per cell, and is the ex¬ 
pected value for Pback when all steps are equiprobable. 
This is because the mean size of an entrance throat is 
larger than the mean size of all the throats. The en¬ 
hanced backward step probability can be regarded as 
a correlation between successively visited throats. As r 
increases, the total available throat area decreases and 
Pback increases, as can be seen in the right panel of fig. 

For very large particles, Pback dominates the walk. 
This can be seen as the ‘lowest order’ effect of connec¬ 
tivity, which we next try to capture within CTRW. 

We examine two methods to modify the CTRW 
model, both introducing a backward bias. In the first, 
we add to P/ and Pt a third distribution, Pg, of step di¬ 
rection relative to the previous step. In this method, the 
particle ‘remembers’ the direction of the previous step, 
and a relative angle is chosen from the non-uniform dis¬ 
tribution Pe{0), e.g. the one in fig. Pe can be calcu¬ 
lated directly from the porous structure for any particle 
radius r, similarly to Pi and Pt (see section |^. A ran¬ 
dom step is then made at the angle 9 relative to the pre¬ 
vious step direction. The waiting time and step length 
at each step are chosen as before. This method intro¬ 
duces a correlation between consecutive steps, which is 
present unavoidably in the diffusion of large particles in 



Fig. 6: Left: the PDF of cos(6), with 9 the angle between 
successive steps of a diffusion process in the porous sample 
for r = ro/100. The singular point at cos{9) = —1 marks the 
finite probability of a backward step Pback = 0.087±0.001. In 
addition, P(—1 < cos(6) < —0.88) = 0. Right: the dependence 
of Pback on r. For large particles, backward steps dominate 
the walk. 


the sample. Testing this method by a set of simulations, 
we find that the universality class transition occurs at 
^(anisotropic) _ ^ ^ Q Q 3 _ ^Joser to the one measured 

in the porous sample. 

The second method is more drastic: at every step 
of the CTRW we construct a new cell, according to the 
structural statistics of the porous sample. We choose 
the cell’s volume, the number of throats, and the throats’ 
areas from the corresponding distributions, derived from 
the sample. Using these and the particle size, we calcu¬ 
late the waiting time for the new cell. Then we choose 
randomly an accessible exit throat. The probability to 
exit through the throat is proportional to its effective 
area. The particle then makes a step in the direction of 
the exit throat. Another cell is then constructed around 
it. The step length is the sum of the radii of these two 
cells. This process is then iterated as many times as 
required. 

A key feature of this method is that the particle 
‘remembers’ the exit direction, the last used cell and 
the exit throat. The latter is used as one of the throats 
in constructing the next cell. If this throat is chosen 
again then the last step is retraced. We refer to this 
method as DA for the two variables that the particle 
‘remembers’ - step direction and throat area. 

Within the DA process, the particle may pass back 
and forth several times through a large throat before 
it moves on and loses memory of this throat. This pro¬ 
cess correlates not only successive steps, as the previous 
method does, but also successive backward steps. 

Using the DA model, we obtain a universality-class 
transition at = 2.50±0.03, in excellent agreement 

with the original simulation results. Thus, this model 
captures much better the particle size-driven topologi¬ 
cal change. 

An important feature of the DA model is that Pback 
is higher than in the porous sample. To understand this. 
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consider two cells in the porous sample, connected by a 
large throat, and connected to the rest of the structure 
by smaller throats. Once the particle enters one of these 
cells, it is likely to move back and forth several times be¬ 
fore it emerges. However, the smaller the throats lead¬ 
ing to this pair of cells, the less likely is the particle 
to enter in the first place. This means that the occur¬ 
rence frequency of such sub-structures, which increase 
^back) is suppressed. In contrast, once a particle passes 
through a relatively large throat in the DA model, then, 
other than this throat, an entire new cell is generated 
for each step. This results in a higher probability that 
the particle oscillates across such a throat. This feature 
appears to compensate for other, more complex, miss¬ 
ing topological features, making the DA a better model 
for the diffusion process. 

As a further investigation of the proposed methods, 
we present the step-length correlation function for the 
different models, all using r = 1.2ro (fig. [^. The corre¬ 
lation in the sample is mainly due to the fact that each 
two consecutive steps enter and exit a certain cell, c. 
If c is small, then the two steps will tend to be short, 
and vice versa, leading to positive correlation. In addi¬ 
tion, a high Pback means that many consecutive steps 
are of exactly the same length. This further increases 
the correlation for larger particles. As expected, the 
traditional CTRW exhibits no correlations. Our first 
adjustment to CTRW, introducing the possibility of a 
backward step, adds correlation. However, since this is 
only a one-step correlation it decays exponentially. The 
increase in correlation within the DA process is because 
of the increased probability for a long sequence of back¬ 
ward steps, as discussed above. As a result, the DA cor¬ 
relation function does not decay exponentially, agreeing 
better with the diffusion simulation in the sample. The 
correlation of the DA at a step distance of one is higher 
than that in the simulated diffusion process because its 
^^back is higher, yet the DA correlation decays faster 
with the number of steps since it lacks the more in¬ 
volved topological correlations. 

6 Conclusions 

To conclude, we compared between two numerical mod¬ 
els of diffusion of a finite size particle in porous media: a 
direct simulation of the diffusion process in a computer¬ 
generated sample, and what is commonly believed to be 
an equivalent CTRW. We first presented a method to 
construct the representative PDFs of the step-length 
and waiting-time. Pi and Pj, given the particle size and 
statistical information about the structure of the porous 
media alone. Using the same particle size dependent 
Pi and Pt in both models, we analysed the transition 



Fig. 7: Step-length correlation function for the four types of 
simulations discussed in the paper, all with r = 1.2ro. Solid 
lines are decaying exponentials. The correlation function of 
the first adjustment to CTRW decays exponentially, and the 
second adjustment (DA) - more slowly. 

from regular to anomalous diffusion. We showed that 
the the two models give different results - while the 
CTRW simulations follow the theoretical prediction, up 
to hnite-time effects, with a transition to sub-diffusion 
at « 2, the diffusion simulations in the confined ge¬ 
ometry exhibit a transition at (3 ~ 2.5. We established 
that the difference stems from changes to the effective 
connectivity available to the particle with increasing 
size. This particle size-driven change in connectivity is 
not taken into consideration in the CTRW model. We 
supported this conclusion by investigation of the time 
average MSD and by showing that the transition in 
universality class occurs at the same range of parti¬ 
cle sizes that corresponds to the percolation transition. 
Our findings show unequivocally that the discrepancy 
between the two models is not due to different waiting¬ 
time distributions, since using identical distributions do 
lead to different universality classes. 

It is important to comment on the range of valid¬ 
ity of our results. Increasing the particle size can be re¬ 
garded, alternatively, as shrinking the porous structure, 
while keeping the particle size unchanged. Evidently, 
the smaller the pore space the more restricted the dif¬ 
fusing particle is and the larger the discrepancy between 
the simulated diffusion and the CTRW. Thus, by start¬ 
ing from a medium with a very large pore space, we 
established that our results hold true for a wide range 
of porous media with lower porosity. 

Wong et al. m studied experimentally a related 
process of trace particles diffusing in biological net¬ 
works of entangled F-actin filaments. There, the dif¬ 
fusion of particles, of size comparable to the typical 
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network mesh size, is sub-linear and Pt decays alge¬ 
braically. The universality class they observe is a func¬ 
tion of the particle-to-mesh size ratio, in agreement 
with our results. However, they do not observe size- 
driven topological effects and their values of a and /3 
are in good agreement with the CTRW model. This is 
because of the flexibility of the gel-like network, which 
allows trapped particles to eventually escape by deform¬ 
ing the filament network, a phenomenon also modelled 
recently in [25] . The rigidity of the structure considered 
here precludes this particle escape mechanism and is the 
reason for this difference. A potentially related process 
of large particles, diffusing in rigid porous media, was 
studied experimentally in |50|. That study focused on 
hydrodynamic in-pore effects, which might be interest¬ 
ing to eventually include in our model. 

The main advantages of the CTRW model are that 
it overcomes potential hnite-size problems and is less 
demanding computationally. However, as we have demon¬ 
strated here, this is achieved at the expense of ignor¬ 
ing topological information about local connectivity. To 
preserve the advantages of CTRW, these need to be 
taken into consideration. To this end, we introduced two 
anisotropic CTRW models. One includes memory of the 
last step direction and a non-uniform distribution of 
step direction. The other, the DA model, adds memory 
of the area of the last throat visited, effectively corre¬ 
lating successive backward steps. The DA model shows 
a universality-class transition at = 2.50±0.03, in 

good agreement with the one measured in our simula¬ 
tions of the diffusion process in the conhned geometry 
of a porous medium. We conclude that the DA model is 
a better alternative to the traditional CTRW for mod¬ 
elling diffusion of hnite size particles in such media. It 
combines the CTRW advantages, overcoming the finite¬ 
ness of the sample and convenience of application, with 
a better capturing of topological effects. 
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